#include "mat.h"


void mat::errorifpointerisnull(void)
{
    if (rawmatptr == NULL)
    {
        std::cout << "Error in 'mat' object: cannot perform the operation (matrix is undefined)" << std::endl;
        abort();
    }
}

void mat::errorifinvalidated(void)
{
    if (rawmatptr != NULL && rawmatptr->getDofManager()->isManaged() && rawmatptr->getmeshnumber() != universe::mymesh->getmeshnumber())
    {
        std::cout << "Error in 'mat' object: matrix cannot be used anymore (invalidated by hp-adaptivity)" << std::endl;
        abort();
    }
}

mat::mat(long long int matsize, intDenseMatrix rowadresses, intDenseMatrix coladresses, denseMatrix vals)
{
    rawmatptr = std::shared_ptr<rawmat>(new rawmat(std::shared_ptr<DofManager>(new DofManager(matsize))));
    rawmatptr->accumulate(rowadresses, coladresses, vals);
    rawmatptr->process();
    rawmatptr->clearfragments();
}

mat::mat(formulation myformulation, intDenseMatrix rowadresses, intDenseMatrix coladresses, denseMatrix vals)
{
    rawmatptr = std::shared_ptr<rawmat>(new rawmat(myformulation.getDofManager()));
    rawmatptr->accumulate(rowadresses, coladresses, vals);
    rawmatptr->process();
    rawmatptr->clearfragments();
}

long long int mat::countRows(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->countRows(); }
long long int mat::countColumns(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->countColumns(); }
        
long long int mat::countNNZ(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->countNNZ(); }

void mat::permute(intDenseMatrix rowpermute, intDenseMatrix colpermute)
{
    if (rowpermute.count() != countRows() || colpermute.count() != countColumns())
    {
        std::cout << "Error in 'mat' object: unexpected argument size for permutation" << std::endl;
        abort();
    }

    Mat permutedmat;
    
    IS rowpermutis, colpermutis;
    ISCreateGeneral(PETSC_COMM_SELF, rowpermute.count(), rowpermute.getValues(), PETSC_USE_POINTER, &rowpermutis);
    ISCreateGeneral(PETSC_COMM_SELF, colpermute.count(), colpermute.getValues(), PETSC_USE_POINTER, &colpermutis);
    
    MatPermute(getPetsc(), rowpermutis, colpermutis, &permutedmat);
    
    rawmatptr = std::shared_ptr<rawmat>(new rawmat(rawmatptr->getDofManager(), permutedmat));
}

void mat::removeConstraints(void) { errorifpointerisnull(); errorifinvalidated(); rawmatptr->removeConstraints(); };

void mat::reuseFactorization(void) { errorifpointerisnull(); errorifinvalidated(); rawmatptr->reuselu(); }

std::shared_ptr<rawmat> mat::getPointer(void)
{
    errorifinvalidated();
    return rawmatptr;
}
        
Mat mat::getPetsc(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getPetsc(); }   

void mat::print(void) { errorifpointerisnull(); errorifinvalidated(); rawmatptr->print(); }

mat mat::copy(void)
{
    Mat A = getPetsc();
    Mat output;
    MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &output);
    return mat(std::shared_ptr<rawmat>(new rawmat(  rawmatptr->getDofManager(), output  )));
}


mat mat::operator+(void) { return copy(); }
mat mat::operator-(void) { return *this*-1; }

mat mat::operator*(double input)
{
    Mat A = getPetsc();
    Mat output;
    MatDuplicate(A, MAT_SHARE_NONZERO_PATTERN, &output);
    MatAXPY(output, input, A, SAME_NONZERO_PATTERN);
    return mat(std::shared_ptr<rawmat>(new rawmat(  rawmatptr->getDofManager(), output  )));
}

mat mat::operator/(double input) { return *this*(1.0/input); }

mat mat::operator*(mat input)
{
    Mat A = getPetsc();
    Mat B = input.getPetsc();
    Mat output;
    MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &output);
    return mat(std::shared_ptr<rawmat>(new rawmat(  rawmatptr->getDofManager(), output  )));
}

mat mat::operator+(mat input)
{
    mat copied = copy();
    Mat Y = copied.getPetsc();
    Mat X = input.getPetsc();
    MatAXPY(Y, 1, X, DIFFERENT_NONZERO_PATTERN);
    return copied;
}

mat mat::operator-(mat input)
{
    mat copied = copy();
    Mat Y = copied.getPetsc();
    Mat X = input.getPetsc();
    MatAXPY(Y, -1, X, DIFFERENT_NONZERO_PATTERN);
    return copied;
}

vec mat::operator*(vec input)
{
    Mat A = getPetsc();
    Vec x = input.getPetsc();
    Vec output;
    VecDuplicate(x, &output);
    MatMult(A, x, output);
    return vec(std::shared_ptr<rawvec>(new rawvec(  rawmatptr->getDofManager(), output  )));
}


mat operator*(double inputdouble, mat inputmat) { return inputmat*inputdouble; }



